home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / interpolation / akima.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  9.4 KB  |  390 lines
  1. /* interpolation/akima.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <math.h>
  25. #include <gsl/gsl_errno.h>
  26. #include "integ_eval.h"
  27. #include <gsl/gsl_interp.h>
  28.  
  29. typedef struct
  30. {
  31.   double * b;
  32.   double * c;
  33.   double * d;
  34.   double * _m;
  35. } akima_state_t;
  36.  
  37.  
  38. /* common creation */
  39. static void *
  40. akima_alloc (size_t size)
  41. {
  42.   akima_state_t *state = (akima_state_t *) malloc (sizeof (akima_state_t));
  43.   
  44.   if (state == NULL)
  45.     {
  46.       GSL_ERROR_NULL("failed to allocate space for state", GSL_ENOMEM);
  47.     }
  48.   
  49.   state->b = (double *) malloc (size * sizeof (double));
  50.   
  51.   if (state->b == NULL)
  52.     {
  53.       free (state);
  54.       GSL_ERROR_NULL("failed to allocate space for b", GSL_ENOMEM);
  55.     }
  56.   
  57.   state->c = (double *) malloc (size * sizeof (double));
  58.   
  59.   if (state->c == NULL)
  60.     {
  61.       free (state->b);
  62.       free (state);
  63.       GSL_ERROR_NULL("failed to allocate space for c", GSL_ENOMEM);
  64.     }
  65.   
  66.   state->d = (double *) malloc (size * sizeof (double));
  67.   
  68.   if (state->d == NULL)
  69.     {
  70.       free (state->c);
  71.       free (state->b);
  72.       free (state);
  73.       GSL_ERROR_NULL("failed to allocate space for d", GSL_ENOMEM);
  74.     }
  75.  
  76.   state->_m = (double *) malloc ((size + 4) * sizeof (double));
  77.  
  78.   if (state->_m == NULL)
  79.     {
  80.       free (state->d);
  81.       free (state->c);
  82.       free (state->b);
  83.       free (state);
  84.       GSL_ERROR_NULL("failed to allocate space for _m", GSL_ENOMEM);
  85.     }
  86.   
  87.   return state;
  88. }
  89.  
  90.  
  91. /* common calculation */
  92. static void
  93. akima_calc (const double x_array[], double b[],  double c[],  double d[], size_t size, double m[])
  94. {
  95.   size_t i;
  96.  
  97.   for (i = 0; i < (size - 1); i++)
  98.     {
  99.       const double NE = fabs (m[i + 1] - m[i]) + fabs (m[i - 1] - m[i - 2]);
  100.       if (NE == 0.0)
  101.     {
  102.       b[i] = m[i];
  103.       c[i] = 0.0;
  104.       d[i] = 0.0;
  105.     }
  106.       else
  107.     {
  108.       const double h_i = x_array[i + 1] - x_array[i];
  109.       const double NE_next = fabs (m[i + 2] - m[i + 1]) + fabs (m[i] - m[i - 1]);
  110.       const double alpha_i = fabs (m[i - 1] - m[i - 2]) / NE;
  111.       double alpha_ip1;
  112.       double tL_ip1;
  113.       if (NE_next == 0.0)
  114.         {
  115.           tL_ip1 = m[i];
  116.         }
  117.       else
  118.         {
  119.           alpha_ip1 = fabs (m[i] - m[i - 1]) / NE_next;
  120.           tL_ip1 = (1.0 - alpha_ip1) * m[i] + alpha_ip1 * m[i + 1];
  121.         }
  122.       b[i] = (1.0 - alpha_i) * m[i - 1] + alpha_i * m[i];
  123.       c[i] = (3.0 * m[i] - 2.0 * b[i] - tL_ip1) / h_i;
  124.       d[i] = (b[i] + tL_ip1 - 2.0 * m[i]) / (h_i * h_i);
  125.     }
  126.     }
  127. }
  128.  
  129.  
  130. static int
  131. akima_init (void * vstate, const double x_array[], const double y_array[],
  132.             size_t size)
  133. {
  134.   akima_state_t *state = (akima_state_t *) vstate;
  135.  
  136.   double * m = state->_m + 2; /* offset so we can address the -1,-2
  137.                                  components */
  138.  
  139.   size_t i;
  140.   for (i = 0; i <= size - 2; i++)
  141.     {
  142.       m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
  143.     }
  144.   
  145.   /* non-periodic boundary conditions */
  146.   m[-2] = 3.0 * m[0] - 2.0 * m[1];
  147.   m[-1] = 2.0 * m[0] - m[1];
  148.   m[size - 1] = 2.0 * m[size - 2] - m[size - 3];
  149.   m[size] = 3.0 * m[size - 2] - 2.0 * m[size - 3];
  150.   
  151.   akima_calc (x_array, state->b, state->c, state->d, size, m);
  152.   
  153.   return GSL_SUCCESS;
  154. }
  155.  
  156.  
  157. static int
  158. akima_init_periodic (void * vstate,
  159.                      const double x_array[],
  160.                      const double y_array[],
  161.                      size_t size)
  162. {
  163.   akima_state_t *state = (akima_state_t *) vstate;
  164.   
  165.   double * m = state->_m + 2; /* offset so we can address the -1,-2
  166.                                  components */
  167.  
  168.   size_t i;
  169.   for (i = 0; i <= size - 2; i++)
  170.     {
  171.       m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
  172.     }
  173.   
  174.   /* periodic boundary conditions */
  175.   m[-2] = m[size - 1 - 2];
  176.   m[-1] = m[size - 1 - 1];
  177.   m[size - 1] = m[0];
  178.   m[size] = m[1];
  179.   
  180.   akima_calc (x_array, state->b, state->c, state->d, size, m);
  181.  
  182.   return GSL_SUCCESS;
  183. }
  184.  
  185. static void
  186. akima_free (void * vstate)
  187. {
  188.   akima_state_t *state = (akima_state_t *) vstate;
  189.  
  190.   free (state->b);
  191.   free (state->c);
  192.   free (state->d);
  193.   free (state->_m);
  194.   free (state);
  195. }
  196.  
  197.  
  198. static
  199. int
  200. akima_eval (const void * vstate,
  201.             const double x_array[], const double y_array[], size_t size,
  202.             double x,
  203.             gsl_interp_accel * a,
  204.             double *y)
  205. {
  206.   const akima_state_t *state = (const akima_state_t *) vstate;
  207.  
  208.   size_t index;
  209.   
  210.   if (a != 0)
  211.     {
  212.       index = gsl_interp_accel_find (a, x_array, size, x);
  213.     }
  214.   else
  215.     {
  216.       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
  217.     }
  218.   
  219.   /* evaluate */
  220.   {
  221.     const double x_lo = x_array[index];
  222.     const double delx = x - x_lo;
  223.     const double b = state->b[index];
  224.     const double c = state->c[index];
  225.     const double d = state->d[index];
  226.     *y = y_array[index] + delx * (b + delx * (c + d * delx));
  227.     return GSL_SUCCESS;
  228.   }
  229. }
  230.  
  231.  
  232. static int
  233. akima_eval_deriv (const void * vstate,
  234.                   const double x_array[], const double y_array[], size_t size,
  235.                   double x,
  236.                   gsl_interp_accel * a,
  237.                   double *dydx)
  238. {
  239.   const akima_state_t *state = (const akima_state_t *) vstate;
  240.  
  241.   size_t index;
  242.  
  243.   DISCARD_POINTER(y_array); /* prevent warning about unused parameter */
  244.   
  245.   if (a != 0)
  246.     {
  247.       index = gsl_interp_accel_find (a, x_array, size, x);
  248.     }
  249.   else
  250.     {
  251.       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
  252.     }
  253.   
  254.   /* evaluate */
  255.   {
  256.     double x_lo = x_array[index];
  257.     double delx = x - x_lo;
  258.     double b = state->b[index];
  259.     double c = state->c[index];
  260.     double d = state->d[index];
  261.     *dydx = b + delx * (2.0 * c + 3.0 * d * delx);
  262.     return GSL_SUCCESS;
  263.   }
  264. }
  265.  
  266.  
  267. static
  268. int
  269. akima_eval_deriv2 (const void * vstate,
  270.                    const double x_array[], const double y_array[], size_t size,
  271.                    double x,
  272.                    gsl_interp_accel * a,
  273.                    double *y_pp)
  274. {
  275.   const akima_state_t *state = (const akima_state_t *) vstate;
  276.  
  277.   size_t index;
  278.  
  279.   DISCARD_POINTER(y_array); /* prevent warning about unused parameter */
  280.  
  281.   if (a != 0)
  282.     {
  283.       index = gsl_interp_accel_find (a, x_array, size, x);
  284.     }
  285.   else
  286.     {
  287.       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
  288.     }
  289.   
  290.   /* evaluate */
  291.   {
  292.     const double x_lo = x_array[index];
  293.     const double delx = x - x_lo;
  294.     const double c = state->c[index];
  295.     const double d = state->d[index];
  296.     *y_pp = 2.0 * c + 6.0 * d * delx;
  297.     return GSL_SUCCESS;
  298.   }
  299. }
  300.  
  301.  
  302. static
  303. int
  304. akima_eval_integ (const void * vstate,
  305.                   const double x_array[], const double y_array[], size_t size,
  306.                   gsl_interp_accel * acc,
  307.                   double a, double b,
  308.                   double * result)
  309. {
  310.   const akima_state_t *state = (const akima_state_t *) vstate;
  311.  
  312.   size_t i, index_a, index_b;
  313.  
  314.   if (acc != 0)
  315.     {
  316.       index_a = gsl_interp_accel_find (acc, x_array, size, a);
  317.       index_b = gsl_interp_accel_find (acc, x_array, size, b);
  318.     }
  319.   else
  320.     {
  321.       index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
  322.       index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
  323.     }
  324.   
  325.   *result = 0.0;
  326.  
  327.   /* interior intervals */
  328.   
  329.   for(i=index_a; i<=index_b; i++) {
  330.     const double x_hi = x_array[i + 1];
  331.     const double x_lo = x_array[i];
  332.     const double y_lo = y_array[i];
  333.     const double dx = x_hi - x_lo;
  334.     if(dx != 0.0) {
  335.  
  336.       if (i == index_a || i == index_b)
  337.         {
  338.           double x1 = (i == index_a) ? a : x_lo;
  339.           double x2 = (i == index_b) ? b : x_hi;
  340.           *result += integ_eval (y_lo, state->b[i], state->c[i], state->d[i],
  341.                                  x_lo, x1, x2);
  342.         }
  343.       else
  344.         {
  345.           *result += dx * (y_lo 
  346.                            + dx*(0.5*state->b[i] 
  347.                                  + dx*(state->c[i]/3.0 
  348.                                        + 0.25*state->d[i]*dx)));
  349.         }
  350.     }
  351.     else {
  352.       *result = 0.0;
  353.       return GSL_FAILURE;
  354.     }
  355.   }
  356.   
  357.   return GSL_SUCCESS;
  358. }
  359.  
  360.  
  361. static const gsl_interp_type akima_type = 
  362. {
  363.   "akima", 
  364.   5,
  365.   &akima_alloc,
  366.   &akima_init,
  367.   &akima_eval,
  368.   &akima_eval_deriv,
  369.   &akima_eval_deriv2,
  370.   &akima_eval_integ,
  371.   &akima_free
  372. };
  373.  
  374. const gsl_interp_type * gsl_interp_akima = &akima_type;
  375.  
  376. static const gsl_interp_type akima_periodic_type = 
  377. {
  378.   "akima-periodic", 
  379.   5,
  380.   &akima_alloc,
  381.   &akima_init_periodic,
  382.   &akima_eval,
  383.   &akima_eval_deriv,
  384.   &akima_eval_deriv2,
  385.   &akima_eval_integ,
  386.   &akima_free
  387. };
  388.  
  389. const gsl_interp_type * gsl_interp_akima_periodic = &akima_periodic_type;
  390.